home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / Helix-turn-helix ƒ / THINK C source code / hth.c next >
Encoding:
C/C++ Source or Header  |  1993-12-28  |  11.2 KB  |  444 lines  |  [TEXT/KAHL]
  1. /*
  2.     hth.c
  3.     v 1.0.5
  4.     28 December 1993
  5.     
  6.     This simple program predicts whether a protein contains a helix-turn-
  7.         helix motif, using the method of:
  8.     
  9.     Dodd, I. B., and J. B. Egan. 1990.  Improved detection of helix-turn-
  10.         helix DNA-binding motifs in protein sequences. Nucleic Acids Res. 
  11.         18:5019-5026.
  12.  
  13.     This code written and donated to the public domain by:
  14.         Conrad Halling
  15.         Department of Molecular Genetics and Cell Biology
  16.         University of Chicago
  17.         920 E 58th St
  18.         Chicago IL 60637
  19.         
  20.         current address:
  21.         
  22.         Monsanto Co.
  23.         700 Chesterfield Parkway North
  24.         St. Louis MO 63198
  25.         e-mail: chhall@bb1t.monsanto.com
  26.     
  27.     How to compile this program:
  28.  
  29.         This program is written in ANSI C using THINK C 5.0.4.
  30.         On the UNIX system at the University of Chicago (running SunOS Release 4.1.1),
  31.             this code will not compile using cc but will compile using either gcc or acc.
  32.         For example, to compile this program, you would type
  33.             acc -o hth hth.c <return>
  34.         This means, "start the acc (ANSI C Compiler) program, send
  35.             the output to a file called "hth", and take the input
  36.             from the file "hth.c".
  37.         When the program is compiled, you run it by typing "hth"
  38.             at the prompt.
  39.  
  40.     tabs = 4
  41.  
  42.         When using vi under UNIX, you can set the tabs to 4 by opening the
  43.             file, typing escape (to go into command mode), colon (":") to go
  44.             to the command line, and "set tabstop=4" (without the quotes).
  45.     
  46.     Format of input protein sequence:
  47.  
  48.         One protein sequence per file
  49.         Single-letter code in upper case and/or lower case
  50.         White space characters (space, tab, return, etc.) are ignored
  51.         The program will abort if an invalid character is found
  52. */
  53.  
  54. #include <ctype.h>
  55. #include <limits.h>
  56. #include <stddef.h>
  57. #include <stdlib.h>
  58. #include <stdio.h>
  59. #include <string.h>
  60.  
  61. #ifndef __HTH__
  62.     #define __HTH__
  63.     
  64.     #ifndef TRUE
  65.         #define TRUE                 1
  66.     #endif
  67.     
  68.     #ifndef FALSE
  69.         #define FALSE                 0
  70.     #endif
  71.  
  72.     #define AMINO_ACIDS_COUNT        20
  73.     #define MAX_SEQUENCE_LENGTH     20000
  74.     
  75.     #define WINDOW_SIZE                22        /* These values from Table 3    */
  76.     #define NON_HTH_MEAN_SCORE       238.71    /* of Dodd and Egan (1991)        */
  77.     #define NON_HTH_STD_DEV           293.61
  78.  
  79.     /*
  80.         Errors
  81.     */
  82.     
  83.     #define NO_ERROR                 0
  84.     #define SEQUENCE_TOO_SHORT         1
  85.     #define OUT_OF_MEMORY             2
  86.     #define QUIT                     3
  87.     #define INVALID_CHAR             4
  88.     
  89.     /*
  90.         Function prototypes
  91.     */
  92.     
  93.     void    DisplayError(
  94.                 short            error );
  95.     void    DisplayResults(
  96.                 double            convertedScore,
  97.                 size_t            maxScorePosition,
  98.                 const char        *sequence );
  99.     short    GetAminoAcid(
  100.                 char            residue );
  101.  
  102. #endif
  103.  
  104.  
  105. const short weightMatrix[ AMINO_ACIDS_COUNT ][ WINDOW_SIZE ] =
  106.     {
  107.     /* A (alanine) */
  108.         -125, -194,  -84,   70,   36,   54,  238,  -15,   77,   26, -194,
  109.         -194,  -56,  -84,   14,   77,  -56,  -56,  -56,   46, -195,   36,
  110.     /* C (cysteine) */
  111.          -64,  -64,  -63,  -63,  -64,  -64,  -64,  -64,  -64,   47,   47,
  112.          -63,  -63,  -64,  -64,  -64,  -64,  -64,  -64,  -63,  -64,   47,
  113.     /* D (aspartate) */
  114.         -156, -154, -156, -154,  109, -156, -156,  109, -154, -156,    6,
  115.         -156, -154,  -85, -156, -156, -156, -156, -154, -154, -156,  -85,
  116.     /* E (glutamate) */
  117.          -31,   -9, -171,   70,  156, -171, -171,  107,   50,  -60,  -60,
  118.         -171,  -60,   78,   86, -171, -171, -101, -171, -170,   86,    9,
  119.     /* F (phenylalanine) */
  120.           10, -130,   10, -130, -130,   10, -130, -129, -130,  102, -130,
  121.         -130, -130, -130, -129, -130, -130, -129, -129, -129,  180, -130,
  122.     /* G (glycine) */
  123.           30,    5, -190,  -51, -191, -191,   18, -191, -191, -191,  202,
  124.         -191,  -10, -191,    5, -190, -191,  -80, -190, -190, -191,  -51,
  125.     /* H (histidine ) */
  126.           62,   33,  -76,  -76,   -7,  -78,  -78,   33,   -7,  -78,   84,
  127.          -78,   33,   33,  -78,   -7,  -78,   -7,   62,   84,  -78,   -7,
  128.     /* I (isoleucine ) */
  129.           75, -156,  101,  -45,  -86,  116, -156,  -16,   65,  -16, -156,
  130.          128, -156,  -86, -156, -155,  188, -155,  -16,   53,  122, -155,
  131.     /* K (lysine ) */
  132.          -31,  -31,   10,   70,   79, -170, -170,   94,   70, -171,   -9,
  133.         -100, -100,   25, -100, -170, -171,   -9,   38,  -31,   -9,  101,
  134.     /* L (leucine) */
  135.           66, -212,   72, -213, -212,  144, -213, -102,   37,  132, -213,
  136.           97, -213, -142, -212, -212,   97, -212, -212,   37,   88, -213,
  137.     /* M (methionine) */
  138.          122,  -74,   -3,  -73,  -73,  -73,  -74,  -74,   88,  122,  -73,
  139.          158,  -74,  -74,   -3,  -74,   -3,  -74,  -74,   -3,  -74,  -73,
  140.     /* N (asparagine) */
  141.         -137,   72, -137, -136, -137, -137, -137,  -67, -136, -136,  128,
  142.         -137,   72, -136,    2,  -67, -137,    2,   84, -137, -137,  104,
  143.     /* P (proline) */
  144.         -156,   23, -157, -156, -157, -157, -157, -157, -157, -157, -157,
  145.         -157,  -46,  101,   39, -157, -157, -157, -157, -157, -157,  -46,
  146.     /* Q (glutamine) */
  147.          -60, -130,  175,   90,  110, -131, -131,   90,   78, -131, -131,
  148.          -60, -130,  154,   65,  119, -131,  -20,  119,   31,  -20,   90,
  149.     /* R (arginine) */
  150.           65,   76,  110,   65,    7, -155, -154,  123,   76, -155, -154,
  151.         -155, -154,  129,   54,   40, -155,  129,  179,  -45, -155,  123,
  152.     /* S (serine) */
  153.         -118,   96, -188,   21,  -48, -187,   -8, -187, -118, -118,  -77,
  154.         -188,  174, -187,  135,  -26, -188,  150,  -77, -188, -187,  -26,
  155.     /* T (threonine) */
  156.           11,  149,  -59,   80, -169,   -8, -170,  -99,  -99,  -30, -170,
  157.           -8,  131,  -30,  -59,  198, -170,  -30, -169, -170,  -30,  -59,
  158.     /* V (valine) */
  159.           17,  -67, -177, -177, -108,  100, -178, -178, -108,   71, -178,
  160.          160, -178,  -16,  -67,  -67,  169, -178,   17,   31,   17, -178,
  161.     /* W (tryptophan) */
  162.           44,  -26,  -26,  -26,  -26,  -26,  -26,  -26,  -26,  -26,  -26,
  163.          -26,  -26,  -25,  -26,  -25,  -26,  -26,   44,  279,  -26,  -26,
  164.     /* Y (tyrosine) */
  165.          -40, -110,   30,    1, -110, -109, -110, -110,  -40,   30, -110,
  166.         -109, -109,  -40, -110,  -40, -110,  162,   52,   86, -110, -110
  167.     };
  168.  
  169. const char aminoAcidsString[] = "ACDEFGHIKLMNPQRSTVWYacdefghiklmnpqrstvwy";
  170.  
  171.  
  172. int main( void )
  173.     {
  174.     char        fileName[ FILENAME_MAX ],
  175.                 format[ 12 ],
  176.                 residue,
  177.                 *sequence;
  178.     int            resultsDisplayed = FALSE,
  179.                 theChar;
  180.     short        aminoAcid,
  181.                 fileNameEntered,
  182.                 maxScore,
  183.                 status,
  184.                 tempScore;
  185.     size_t        i,
  186.                 j,
  187.                 length,
  188.                 maxWindowPosition,
  189.                 maxScorePosition,
  190.                 position,
  191.                 sequenceLength;
  192.     double        convertedScore;
  193.     FILE        *sequenceFile;
  194.     
  195.     status = NO_ERROR;
  196.     
  197.     printf( "Welcome to hth, a program that predicts whether a protein contains\n" );
  198.     printf( "a helix-turn-helix motif.\n\n" );
  199.     printf( "For more information, please read\n\n" );
  200.     printf( "    Dodd, I. B., and J. B. Egan. 1990.  Improved detection of\n" );
  201.     printf( "        helix-turn-helix DNA-binding motifs in protein sequences.\n" );
  202.     printf( "        Nucleic Acids Res. 18:5019-5026.\n\n" );
  203.     printf( "Please be sure that your protein sequence is in \"plain\" format\n" );
  204.     printf( "using the single-letter code.\n\n" );
  205.  
  206.     /*
  207.         Allocate memory for the protein sequence.
  208.     */
  209.  
  210.     sequence = ( char * ) malloc( MAX_SEQUENCE_LENGTH * sizeof( char ) );
  211.     if ( NULL == sequence )
  212.         status = OUT_OF_MEMORY;
  213.  
  214.     while ( NO_ERROR == status )
  215.         {
  216.     
  217.         /*
  218.             Get the name of the sequence file.
  219.         */
  220.  
  221.         fileNameEntered = FALSE;
  222.         while ( TRUE )
  223.             {
  224.             printf( "Name of protein file (q/Q = Quit):  " );
  225.             sprintf( format, "%%%us", FILENAME_MAX );
  226.             scanf( format, fileName );
  227.             length = strlen( fileName );
  228.             if ( 0 != length )
  229.                 {
  230.                 /*
  231.                     Check for quit command.
  232.                 */
  233.                 
  234.                 if ( ( 1 == length ) &&
  235.                     ( ( fileName[ 0 ] == 'q' ) || ( fileName[ 0 ] == 'Q' ) ) )
  236.                     {
  237.                     status = QUIT;
  238.                     break;
  239.                     }
  240.                 else
  241.                     {
  242.                     /*
  243.                         Open the sequence file.
  244.                     */
  245.             
  246.                     sequenceFile = fopen( fileName, "r" );
  247.                     if ( NULL == sequenceFile )
  248.                         printf( "File not found!\n\n" );
  249.                     else
  250.                         break;
  251.                     }
  252.                 }
  253.             }
  254.  
  255.         if ( NO_ERROR == status )
  256.             {
  257.             /*
  258.                 Read the sequence from sequenceFile into sequence[].
  259.             */
  260.             
  261.             i = 0;
  262.             while ( i < MAX_SEQUENCE_LENGTH )
  263.                 {
  264.                 /*
  265.                     Stop reading at end of file.
  266.                 */
  267.                 
  268.                 theChar = getc( sequenceFile );
  269.                 if ( feof( sequenceFile ) )
  270.                     break;
  271.     
  272.                 /*
  273.                     Skip white space characters.
  274.                 */
  275.                 
  276.                 if ( isspace( theChar ) )
  277.                     continue;
  278.     
  279.                 /*
  280.                     GetAminoAcid will return AMINO_ACIDS_COUNT if the character
  281.                         is not found in AminoAcidsString[].
  282.                 */
  283.     
  284.                 if ( AMINO_ACIDS_COUNT == GetAminoAcid( ( char ) theChar ) )
  285.                     {
  286.                     status = INVALID_CHAR;
  287.                     break;
  288.                     }
  289.     
  290.                 /*
  291.                     The character is valid; add it to the string.
  292.                 */
  293.     
  294.                 sequence[ i++ ] = ( char ) theChar;
  295.                 }
  296.     
  297.             if ( NO_ERROR == status )
  298.                 sequence[ i++ ] = 0x00;
  299.         
  300.             fclose( sequenceFile );
  301.             }
  302.     
  303.         if ( NO_ERROR == status )
  304.             {
  305.             sequenceLength = strlen( sequence );
  306.             if ( sequenceLength < WINDOW_SIZE )
  307.                 status = SEQUENCE_TOO_SHORT;
  308.             }
  309.         
  310.         if ( NO_ERROR == status )
  311.             {
  312.             /*
  313.                 Calculate the highest score for the sequence.
  314.             */
  315.     
  316.             resultsDisplayed = FALSE;
  317.             maxScore = SHRT_MIN;    /* defined in <limits.h> */
  318.             maxScorePosition = 0;
  319.             maxWindowPosition = sequenceLength - WINDOW_SIZE;
  320.             for ( i = 0; i < maxWindowPosition; i++ )
  321.                 {
  322.                 tempScore = 0;
  323.                 for ( j = 0; j < WINDOW_SIZE; j++ )
  324.                     {
  325.                     position = i + j;
  326.                     residue = sequence[ position ];
  327.                     aminoAcid = GetAminoAcid( residue );
  328.                     tempScore += weightMatrix[ aminoAcid ] [ j ];
  329.                     }
  330.                 if ( tempScore > maxScore )
  331.                     {
  332.                     maxScore = tempScore;
  333.                     maxScorePosition = i;
  334.                     convertedScore =
  335.                         ( ( double ) maxScore - NON_HTH_MEAN_SCORE ) /
  336.                             ( NON_HTH_STD_DEV );
  337.                     if ( convertedScore >= 2.5 )
  338.                         {
  339.                         DisplayResults(
  340.                             convertedScore,
  341.                             maxScorePosition,
  342.                             sequence );
  343.                         resultsDisplayed = TRUE;
  344.                         maxScore = SHRT_MIN;
  345.                         }
  346.                     }
  347.                 }
  348.             if ( !resultsDisplayed )
  349.                 DisplayResults( convertedScore, maxScorePosition, sequence );
  350.             }
  351.         if ( ( QUIT != status ) && ( NO_ERROR != status ) )
  352.             {
  353.             DisplayError( status );
  354.             status = NO_ERROR;
  355.             }
  356.         }
  357.     DisplayError( status );
  358.     }
  359.  
  360.  
  361. void DisplayResults(
  362.     double            convertedScore,
  363.     size_t            maxScorePosition,
  364.     const char        *sequence )
  365.     {
  366.     char            maxScoreString[ WINDOW_SIZE + 1 ];
  367.     short            i,
  368.                     percentage;
  369.     size_t            position;
  370.     
  371.     printf(
  372.         "The score is %0.2f at position %lu.\n", 
  373.         convertedScore,
  374.         maxScorePosition + 1 );
  375.  
  376.     for ( i = 0; i < WINDOW_SIZE; i++ )
  377.         {
  378.         position = maxScorePosition + i;
  379.         maxScoreString[ i ] = sequence[ position ];
  380.         }
  381.     maxScoreString[ i ] = 0x00;
  382.  
  383.     printf(
  384.         "The sequence at this position is %s.\n",
  385.         maxScoreString );
  386.  
  387.     if ( convertedScore >= 4.5 )
  388.         percentage = 100;
  389.     else if ( convertedScore >= 4.0 )
  390.         percentage = 90;
  391.     else if ( convertedScore >= 3.5 )
  392.         percentage = 71;
  393.     else if ( convertedScore >= 3.0 )
  394.         percentage = 50;
  395.     else if ( convertedScore >= 2.5 )
  396.         percentage = 25;
  397.  
  398.     if ( convertedScore < 2.5 )
  399.         printf( "This score is not significant.\n\n" );
  400.     else
  401.         {
  402.         printf(
  403.             "This score suggests an approximately %hd%% probability that ",
  404.             percentage );
  405.         printf( "this protein\ncontains a helix-turn-helix motif.\n\n" );
  406.         }
  407.     }
  408.                 
  409.  
  410. short GetAminoAcid(
  411.     char        residue )
  412.     {
  413.     short        i,
  414.                 limit;
  415.     
  416.     limit = strlen( aminoAcidsString );
  417.     for ( i = 0; i < limit; i++ )
  418.         {
  419.         if ( residue == aminoAcidsString[ i ] )
  420.             break;
  421.         }
  422.     if ( i >= AMINO_ACIDS_COUNT )
  423.         i -= AMINO_ACIDS_COUNT;
  424.  
  425.     return ( i );
  426.     }
  427.  
  428.  
  429. void DisplayError(
  430.     short        error )
  431.     {
  432.     
  433.     char        errorString[5][80] =
  434.         {
  435.         "\n\n",
  436.         "The protein sequence is too short to analyze.\n\n",
  437.         "There is insufficient memory to continue.\n\n",
  438.         "Good-bye!\n\n",
  439.         "There is an invalid character in the protein sequence.\n\n"
  440.         };
  441.  
  442.     printf( errorString[ error ] );
  443.     }
  444.